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Abstract Modern cosmological observations allow us to study in great detail the evo- 
lution and history of the large scale structure hierarchy. The fundamental problem 
of accurate constraints on the cosmological parameters, within a given cosmological 
model, requires precise modelling of the observed structure. In this paper we briefly 
review the current most effective techniques of large scale structure simulations, em- 
phasising both their advantages and shortcomings. Starting with basics of the direct 
N-body simulations appropriate to modelling cold dark matter evolution, we then dis- 
cuss the direct-sum technique GRAPE, particle-mesh (PM) and hybrid methods, com- 
bining the PM and the tree algorithms. Simulations of baryonic matter in the Universe 
often use hydrodynamic codes based on both particle methods that discretise mass, 
and grid-based methods. We briefly describe Eulerian grid methods, and also some 
variants of Lagrangian smoothed particle hydrodynamics (SPH) methods. 
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1 Introduction 

In the hierarchical picture of structure formation, small objects collapse first and then 
merge to form larger and larger structures in a complex manner. This formation pro- 
cess reflects on the intricate structure of galaxy clusters, whose properties depend on 
how the thousands of smaller objects that the cluster accretes are destroyed or sur- 
vive within the cluster gravitational potential. These merging events are the source of 
shocks, turbulence and acceleration of relativistic particles in the intracluster medium, 
which, in turn, lead to a redistribution or amplification of magnetic fields, and to the 
acceleration of cosmic rays. In order to model these processes realistically, we need to 
resort to numerical simulations which are capable of resolving and following correctly 
the highly non-linear dynamics. In this paper, we briefly describe the methods which 
are commonly used to simulate galaxy clusters within a cosmological context. 

Usually, choosing the simulation setup is a compromise between the size of the 
region that one has to simulate to fairly represent the object(s) of interest, and the 
resolution needed to resolve the objects at the required level of detail. Typical sizes of 
the simulated volume are a megaparsec scale for an individual galaxy, tens to hundreds 
of megaparsecs for a galaxy population, and several hundreds of megaparsecs for a 
galaxy cluster population. The mass resolution varies from w I0 5 M up to w I0 10 M Q , 
depending on the object studied, while, nowadays, one can typically reach the resolution 
of a few hundred parsec for individual galaxies and above the kiloparsec scale for 
cosmological boxes. 



2 N-Body (pure gravity) 

Over most of the cosmic time of interest for structure formation, the Universe is dom- 
inated by dark matter. The most favourable model turned out to be the so-called 
cold dark matter (CDM) model. The CDM can be described as a collisionless, non- 
relativistic fluid of particles of mass m, position x and momentum p. In an expand- 
ing background Universe (usually described by a Friedmann-Lemaitre model), with 
a = (I + z)~ being the Universe scale factor, x is the comoving position and the 
phase-space distribution function /(x, p, t) of the dark-matter fluid can be described 
by the collisionless Boltzmann (or Vlasov) equation 

S + - E 5W-mV*^=0 (1) 
at ma z op 

coupled with the Poisson equation 

V 2 <£(x, t) = 4nGa 2 [p(x, t) - p{t)] , (2) 

where $ is the gravitational potential and p(t) is the background density. The proper 
mass density 

p(x,i) = J /(x,p,t)d 3 p (3) 

can be inferred by integrating the distribution function over the momenta p = ma 2 x. 

This set of equations represents a high-dimensional problem. It is therefore usually 
solved by sampling the phase-space density by a finite number N of tracer particles. 
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The solution can be found through the equation of motion of the particles (in comoving 
coordinates) , 

^ = -mV$ (4) 

At v ; 

and 

— - -5- (5) 

M ma 2 ' 1 ' 

Introducing the proper peculiar velocity v = ax these equations can be written as 

dv a V<5 

— + v~ = . (6) 

At a a 

The time derivative of the expansion parameter, a, can be obtained from the Friedmann 
equation 

a = h y/i + noia- 1 - 1) + n A (a 2 - 1), (7) 

where we have assumed the dark energy to be equivalent to a cosmological constant. 
For a m ore detai l ed des cription of the underlying cosmology and related issues, see for 
example [Peebles! |l980h or others. 

There are different approaches: to solve directly the motion of the tracer particles, 
or to solve the Poisson equation. Some of the most common methods will be described 
briefly in the following sections. 



2.1 Direct sum (GRAPE, GPU) 

The most direct way to solve the N-body problem is to sum directly the contributions 
of all the individual particles to the gravitational potential 

*(r) = -GV ^ -. (8) 

" (|r-r,|2+ £ 2)5 

In principle, this sum would represent the exact (Newtonian) potential which gener- 
ates the particles' acceleration. As mentioned before, the particles do not represent 
individual dark matter particles, but should be considered as Monte Carlo realisations 
of the mass distribution, and therefore only collective, statistical properties can be 
considered. In such simulations, close encounters between individual particles are irrel- 
evant to the physical problem under consideration, and the gravitational force between 
two particles is smoothed by introducing the gravitational softening e. This softening 
reduces the spurious two-body relaxation which occurs when the number of particles 
in the simulation is not large enough to represent correctly a collisionless fluid. This 
situation however is unavoidable, because the number of dark matter particles in real 
systems is orders of magnitude larger than the number that can be handled in a nu- 
merical simulation. Typically, e is chosen to be 1/20 — 1/50 of the mean inter-particle 
separation within the simulation. In general, this direct-sum approach is considered to 
be the most accurate technique, and is used for problems where superior precision is 
needed. However this method has the disadvantage of being already quite CPU inten- 
sive for even a moderate number of particles, because the computing time is oc N 2 , 
where N is the total number of particles. 

Rather than searching for other software solutions, an alternative approach to solve 
the 7V 2 -bottleneck of the d irect-sum techn ique is the GRAPE (GRAvity PipE) special- 
purpose hardware (see e.g. llto et al.lll993l and related articles). This hardware is based 
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Fig. 1 Schematic illustration of the iBarnes fc Hutl IU986T) oct-trce in two dimensions. The 
particles are first enclosed in a square (root node). This square is then iteratively subdivided 
into four squares of half the size, until exactly one particle is left in each final square (leaves 
of the tree). In the resulting tree s t ructure , each square can be the progenitor of up to four 
siblings. Taken from ISpringel et al.l ll2001bf) . 



on custom chips that compute the gravitational force with a hardwired Plummer force 
law (Eq. [H} . This hardware device thus solves the gravitational N-body problem with 
a direct summation approach at a computational speed which is considerably higher 
than that of traditional processors. 

For the force computation, the particle coordinates are first loaded onto the GRAPE 
board, then the forces for several positions (depending on the number of individual 
GRAPE chips installed in the system) are computed in parallel. In practice, there are 
some technical complications when using the GRAPE system. One is that the hard- 
ware works internally with special fixed-point formats or with limited floating point 
precision (depending on the version of the GRAPE chips used) for positions, accelera- 
tions and masses. This results in a reduced dynamic range compared to the standard 
IEEE floating point arithmetic. Furthermore, the communication time between the 
host computer and the GRAPE system can be an issue in certain circumstances. How- 
ever, newer versions of the GRAPE chips circumvent this problem, and can also be 
co mbined with the t r ee alg o rithms (which a r e described in detail i n the next section), 



_ .Fukushige et 
(|2000h . 



ie tree algorithms (which arc described in d etail m the next section), 
aD i|l99ll ); iMakinol (|l99ll ); lAthanassoula et all (|l998l ); iKawai et al] 



By contrast, the graphic processing unit (GPU) on modern graphic cards now pro- 
vides an alternative tool for high-performance computing. The original purpose of the 
GPU is to serve as a graphics accelerator for speeding up the image processing, thereby 
allowing one to perform simple instructions on multiple data. It has therefore become 
an active area of research to use the GPUs of the individual members of computer 
clusters. Although very specialised, many of those computational algorithms are also 
needed in computational astrophysics, and therefore the GPU can provide significantly 
more computing power than the host system; thereby providing a high performance 
with typically large memory size and at relatively low cost, which represents a valid 
alternative to special purpose hardware li ke GRAPE. For recent applications to astro- 



alternative to special purpose hardware like (jtiAFE. box recent 
physical problems see lSchive et al. and references therein 
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2.2 Tree 

The primary method of solving the N-body problem is a hierarchical multipole expan- 
sion, commonly called a tree algorithm. This method groups distant particles into larger 
cells, allowing their gravity to be accounted for by means of a single multipole force. 
Instead of requiring TV — 1 partial force evaluations per particle, as needed in a direct- 
summation approach, the gravitational force on a single particle can be computed with 
substantially fewer operations, because distant groups are treated as "macro" particles 
in the sum. In this manner the sum usually reduces to Nlog(N) operations. Note how- 
ever that this scaling is only true for homogeneous particle distributions, whereas the 
scaling for strongly inhomogeneous distributions, as present in evolved cosmological 
structures, can be less efficient. 

In practice, the hierarchical grouping that forms the basis of the multipole expan- 
sion is most common ly obtained by a recursive subdivision of space. In the approach of 
iBarnes fc Hut] l|l986f ). a cubical root node is used to encompass the full mass distribu- 



tion; the cube is repeatedly subdivided into eight daughter nodes of half the side-length 
each, until one ends up with 'leaf nodes containing single particles (see Fig.[T|. Forces 
are then obtained by "walking" the tree. In other words, starting at the root node, a 
decision is made as to whether or not the multipole expansion of the node provides 
an accurate enough partial force. If the answer is 'yes', the multipole force is used and 
the walk along this branch of the tree can be terminated; if the answer is 'no', the 
node is "opened", i.e. its daughter nodes are considered in turn. Clearly, the multipole 
expansion is in general appropriate for nodes that are sufficiently small and distant. 
Most commonly one uses a fixed angle (typically « 0.5 rad) as opening criteria. 

It should be noted that the final result of the tree algorithm will in general only 
represent an approximation to the true force. However, the error can be controlled con- 
veniently by modifying the opening criterion for tree nodes, because a higher accuracy 
is obtained by walking the tree to lower levels. Provided that sufficient computational 
resources are invested, the tree force can then be made arbitrarily close to the well- 
specified correct force. Nevertheless evaluating the gravitational force via a tree leads 
to an inherent asymmetry in the interaction between two particles. It is worth men- 
tioning that there are extensions to the standard tree, the so-called fast multipole 
methods, which avoid these asymmetries, and therefore have be tter con s ervati on of 
momentum. For an N-body application of such a technique see iDehnenl (|2000h and 
references therein. However, these methods compute the forces for all the particles at 
every time step and can not take advantage of using individual time steps for different 
particles. 



2.3 Particle-Mesh methods 

The Particle- Mesh (PM) method treats the force as a field quantity by computing it on 
a mesh. Differential operators, such as the Laplacian, are replaced by finite difference 
approximations. Potentials and forces at particle positions are obtained by interpolation 
on the array of mesh-defined values. Typically, such an algorithm is performed in three 
steps. First, the density on the mesh points is computed by assigning densities to the 
mesh from the particle positions. Second, the density field is transformed to Fourier 
space, where the Poisson equation is solved, and the potential is obtained using Green's 
method. Alternatively, the potential can be determined by solving Poisson's equation 
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iteratively with relaxation methods. In a third step the forces for the individual particles 
are obtained by interpolating the derivatives of the potentials to the particle positions. 
Typically, the amount of mesh cells N used corresponds to the number of particles in 
the simulation, so that when structures form, one can have large numbers of particles 
within individual mesh cells, which immediately illustrates the shortcoming of this 
method; namely its limited resolution. On the other hand, the calculation of the Fourier 
transform via a Fast Fourier Transform (FFT) is extremely fast, as it only needs of order 
N log N operations, which is the advantage of this method. Note that here N denotes 
the number of mesh cells. In this approach the computational costs do not depend on 
the details of the particle distribution. Also this method can not take advantage of 
individual time steps, as the forces are always calculated for all particles at every time 
step. 

There are many schemes to assign the mass density to the mesh. The simplest 
method is the "Nearest-Grid-Point" (NGP). Here, each particle is assigned to the 
closest mesh point, and the density at each mesh point is the total mass assigned 
to the point divided by the cell volume. However, this method is rarely used. One 
of its drawbacks is that it gives forces that are discontinuous. The "Cloud-in-a-Cell" 
(CIC) scheme is a better approximation to the force: it distributes every particle over 
the nearest 8 grid cells, and then weighs them by the overlapping volume, which is 
obtained by assuming the particle to have a cubic shape of the same volume as the mesh 
cells. The CIC method gives continuous forces, but discontinuous first derivatives of 
the forces. A more accurate scheme is the "Triangular-Shaped-Cloud" (TSC) method. 
This scheme has an assignment interpolation fu nction that is piecewise qua dratic. In 
three dimensions it employs 27 mesh points fsee lHocknev fc Eastwoodlll988l ). 

In general, one can define the assignment of the density p m on a grid x m with spac- 
ing 8 from the distribution of particles with masses m t and positions Xj, by smoothing 
the particles over n times the grid spacing (h — nS). Therefore, having defined a 
weighting function 

W{xm - Xi) = J W (^^) S(x - Xi, fc)dx, (9) 

where W(x) is 1 for |x| < 0.5 and otherwise, the density p m on the grid can be 
written as 

Pm = ^m l W(x l - x m ). (10) 

i 

The shape function S(x, h) then defines the different schemes. The aforementioned 
NGP, CIC and TSC schemes are equivalent to the choice of 1, 2 or 3 for n and the 
Dirac 8 function 5(x), W{x/h) and 1— |x/ft| for the shape function S(x, h), respectively. 

In real space, the gravitational potential <& can be written as the convolution of the 
mass density with a suitable Green's function g(x): 

*(x) = y.g(x-x')p(x')c!x'. 
For vacuum boundary conditions, for example, the gravitational potential is 

*(*) = -G / r^dx', (12) 
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with G being the gravitational constant. Therefore the Green's function, <?(x) = — G/|x|, 
represents the solution of the Poisson equation V 2 $(x) = 47rGp(x), recalling that 
V^(|x — x'|) _1 = 47r<5(x — x'). By applying the divergence theorem to the integral form 
of the above equation, it is then easy to see that, in spherical coordinates, 

x y2 dv = i s v g) ds = r £ * (9 = - 4 - (i3) 

Periodic boundary conditions are usually used to simulate an "infinite universe" , how- 
ever zero padding can be applied to deal with vacuum boundary conditions. 

In the PM method, the solution to the Poisson equation is easily found in Fourier 
space, where Eg. 1111 becomes a simple multiplication 

<?(k) = 3(k)p(k). (14) 

Note that g(k) has only to be computed once, at the beginning of the simulation. 

After the calculation of the potential via Fast Fourier Transform (FFT) methods, 
the force field f(x) at the position of the mesh points can be obtained by differentiating 
the potential, f (x) = V#(x). This can be done by a finite-difference representation of 
the gradient. In a second order scheme, the derivative with respect to the x coordinate 
at the mesh positions m — k) can be written as 

Jx) _ &i+l,j,k ~ &i-l,j,k / 1t .s 
Ji,j,k - 2h • 1 J 

A fourth order scheme for the derivative would be written as 

Jx) ■!,;./■ :.,-.<■ , l ®i+2,j,k-®i-2,j,k rifi , 

J «. fc ~ 3 2h + 3 Ah ■ [ ' 

Finally, the forces have to be interpolated back to the particle positions as 

f ( Xi ) = V W(x j; - x TO )fm, (17) 



where it is recommended to use the same weighting scheme as for the density assign- 
ment; this ensures pairwise force symmetry between particles and momentum conser- 
vation. 

The advantage of such PM methods is the speed, because the number of operations 
scales with N + N g log(N g ), where N is the number of particles and N g the number 
of mesh points. However, the disadvantage is that the dynamical range is limited by 
Ng, which is usually limited by the available memory. Therefore, particularly for cos- 
mological simulations, adaptive methods are needed to increase the dynamical range 
and follow the formation of individual objects. 

In the Adaptive Mesh Refinement (AMR) techniques, the Poisson equation on 
the refinement meshes can be treated as a Dirichlet boundary problem for which the 
boundary values are obtained by interpolating the gravitational potential from the 
parent grid. In such algorithms, the boundaries of the refinement meshes can have an 
arbitrary shape; this feature narrows the range of solvers that one can use for partial 
differential equation (PDEs). The Poisson equation on t hese meshes can be solved 
using the relaxation method ( Hocknev fc Eastwood! 1 19881 ; IPress et al.l Il992h , which is 



relatively fast and efficient in dealing with complicated boundaries. In this method the 
Poisson equation 

V 2 <£ = p (18) 
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Fig. 2 A slice through the refinement structure (the base grid is not shown) in a /1CDM simu- 
lation (left panel) and the corresponding slice through the particle distribution (middle panel) 
The a rea enclosed by the square is enlarged in the right panel. Taken from iKravtsov et al 



is rewritten in the form of a diffusion equation, 

= V 2 <£ - p. (19) 

or 

The point of the method is that an initial solution guess # relaxes to an equilibrium 
solution (i.e., solution of the Poisson equation) as r — > oo. The finite-difference form 
of Eq. 2 is: 

"•!:;.' = *u* + % ( E - m liA p^ At - c 20 ) 



\nf)=l 



where the summation is performed over a cell's neighbours. Here, A is the actual spatial 
resolution of the solution (potential) , while At is a fictitious time step (not related to 
the actual time integration of the iV-body system). This finite difference method is 
stable when At < A 2 /6. More det ails can be found in P ress et al. (1992) and also 
IKravtsov et al. | 1 19971 ). Fig. [2] from IKravtsov et~aH l|l997l ). shows an example of the 



mesh constructed to calculate the potential in a cosmological simulation. 



2.4 Hybrids (TreePM/P 3 M) 

Hybrid methods can be constructed as a synthesis of the particle-mesh method and the 
tree a lgorithm. In TreePM methods (|Xull995l ; lBode et al.M2000l ; lBaglall2002l ; lBagla fc Ravi 
l2003h the potential is explicitly split in Fourier space into a long-range and a short- 
range part according to $ k = $ k ™ 8 + <j5^ hort , where 

4 ong = <2> k exp(-k 2 r s 2 ), (21) 

with r s describing the spatial scale of the force-split. The long range potential can be 
computed very efficiently with mesh-based Fourier methods. 

The short-range part of the potential can be solved in real space by noting that for 
r s <<iC L the short-range part of the real-space solution of the Poisson equation is given 
by 
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Fig. 3 Force decomposition and force error of the TreePM scheme in the case when two 
stacked meshes are used. The left panel illustrates the strength of the short-range (dot-dashed), 
intermediate-range (thick solid), and long-range (solid) force as a function of distance in a 
periodic box. The spatial scales of the two splits are marked with vertical dashed lines. The 
right panel shows the error distribution of the PM force. The outer matching region exhibits 
a very similar error characteristic as the inner match of tree- and PM-force. In both cases, for 
separations of order the fine or coarse mesh scale (dotted lines), respectively, force errors of 
up to 1 — 2 % arise, but the r.m.s. force e rror stays well below 1 %, and the mean force tracks 
the correct result accurately. Taken from [Sprinecl (20051). 



^short 



w = - G Ef- fc (i)- (22) 



Here is the distance of any particle i to the point x. Thus the short-range force 
can be computed by the tree algorithm, except that the force law is modified by a 
long-range cut-off factor. 

Such hybrid methods can result in a very substantial improvement of the perfor- 
mance compared to ordinary tree methods. In addition one typically gains accuracy in 
the long-range force, which is now basically exact, and not an approximation as in the 
tree method. Furthermore, if r s is chosen to be slightly larger than the mesh scale, force 
anisotropics, that exist in plain PM methods, can be suppressed to essentially arbitrar- 
ily low levels. A TreePM approach also maintains all the most important advantages 
of the tree algorithm, namely its insensitivity to clustering, its essentially unlimited 
dynamical range, and its precise control of the softening scale of the gravitational 

force. 

Fig. 13 shows how the force matching works in the GADGET-2 code ( Springe! 



2005), where such a hybrid method is further extended to two subsequent stacked 
FFTs combined with the tree algorithm. This extension enables one to increase the 
dynamic range, which, in turn, improves the computational speed in high resolution 
simulations of the evolution of galaxy clusters within a large cosmological volume. 

Although used much earlier (because much easier to implement), the P 3 M method 
can be seen as a special case of the TreePM, where the tree is replaced by the direct 
sum. Note that also in the tree algorithm the nearest forces are calculated by a di- 
rect sum, thus the P 3 M approach formally corresponds to extending the di rect sum of 
the tr ee method to the scale where the PM force computation takes over. ICouchmanl 
( 199l|) presented an improved version of the P 3 M method, by allowing spatially adap- 



tive mesh refinements in regions with high particle density (Adaptive P 3 M or AP 3 M). 
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The improvement in the performance made it very attractive and several cosmological 
simul ations were perform ed with this technique, including the Hubble Volume simula- 
tions (jEvrard et alj|2002h . 



2.5 Time-stepping and integration 

The accuracy obtained when evolving the system depends on the size of the time step 
and on the integrator scheme used. Finding the optimum size of time step is not trivial. 
A very simple criterion often used is 



At = oiy/e/\a\ (23) 

where \a\ is the acceleration obtained at the previous time step, e is a length scale, 
which can typically be associated with the gravitational softening, and a is a tolerance 
parameter. More d etails about different time step criteria can be found for example in 
IPower et all (|2003h and references therein. For the integration of the variables (posi- 



tions and velocities) of the system, we only need to integrate first order equations of 
the form y = /(y), e.g. ordinary differential equations (ODEs) with appropriate initial 
conditions. Note that we can first solve this ODE for the velocity v and then treat 
x = v as an independent ODE, at basically no extra cost. 

One can distinguish implicit and explicit methods for propagating the system from 
step n to step n + 1. Implicit methods usually have better properties, however they 
need to solve the system iteratively, which usually requires inverting a matrix which is 
only sparsely sampled, and has the dimension of the total number of the data points, 
namely grid or particle points. Therefore, N-body simulations mostly adopt explicit 
methods. 

The simplest (but never used) method to perform the integration of an ODE is 
called Euler's method; here the integration is just done by multiplying the derivatives 
with the length of the time step. The explicit form of such a method can be written as 

y«+i = y„ + f(y n )At, (24) 

whereas the implicit version is written as 

y n +i =yn + f(y„+i)^t (25) 

Note that in the latter equation y„+i appears on the left and right side, which makes 
it clear why it is called implicit. Obviously the drawback of the explicit method is that 
it assumes that the derivatives (e.g. the forces) do not change during the time step. 

An improvement to this method can be obtained by using the mean derivative 
during the time step, which can be written with the implicit mid-point rule as 

y n +i = yn + f [0.5(y„ + y n+1 )]At. (26) 

An explicit rule using the forces at the next time step is the so-called predictor-corrector 
method, where one first predicts the variables for the next time step 

y°+i = y« + f(y n )At (27) 

and then uses the forces calculated there to correct this prediction (the so-called cor- 
rector step) as 

yn+l = Yn + 0.5[f(y„)+f(y° +1 )]zli. (28) 
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Fig. 4 The upper two rows shows a Kepler problem of high eccentricity evolved with different 
simple time integration schemes, using an equal time step in all cases. Even though the leap-frog 
(upper left panel) and the second order Runge-Kutta (upper right panel) produce comparable 
errors in a single step, the long term stability of the integration is very different. Even a 
computationally much more expensive fourth order Runge-Kutta scheme (middle row), with 
a smaller error per step, performs dramatically worse than the leap-frog in this problem. The 
lower row shows the same using leap-frog schemes with a variable time step from step to step, 
based on the At oc 1/^/jaf criterion commonly employed in cosmological simulations. As a 
result of the variable time steps, the integration is no longer manifestly time reversible, and 
long term secular errors develop. Interestingly, the error in the KDK (Kick-Drift-Kick) variant 
grows four times more slowly than i n the DKD (Dr ift-Kick-Drift) variant, despite being of 
equal computational cost. Taken from Springcl (2005). 
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This method is accurate to second order. 

In fact, all these methods are special cases of the so-called Runge-Kutta method 
(RK), which achieves the accuracy of a Taylor series approach without requiring the 
calculation of higher order derivatives. The price one has to pay is that the derivatives 
(e.g. forces) have to be calculated at several points, effectively splitting the interval At 
into special subsets. For example, a second order RK scheme can be constructed by 

ki = f (y„) (29) 

k 2 = f(y« +kiAfc) (30) 

y n +i = yn + 0.5(ki + k 2 )At. (31) 

In a fourth order RK scheme, the time interval At also has to be subsampled to calculate 
the mid-points, e.g. 

ki = f(y„,tn) (32) 
k 2 = f (yn + kiAt/2, t n + At/2) (33) 
k 3 = f (y« + k 2 At/2, t n + At/2) (34) 
k i = f(y n + k 3 At/2,t n + At) (35) 

-«=-+(t + t + t + t)^ <*> 

More details on how to constru ct the coefficient for an n-th order RK scheme are given 
e.g. lChapra fc Canalel [|l997h . 
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Another possibility is to use the so-called leap-frog method, where the derivatives 
(e.g. forces) and the positions are shifted in time by half a time step. This feature can 
be used to integrate directly the second order ODE of the form x = f (x). Depending 
on whether one starts with a drift (D) of the system by half a time step or one uses 
the forces at the actual time to propagate the system (kick, K), one obtains a KDK 
version 

v„+i/2 = v„ + f(xn)At/2 (37) 
x n+ i = x n + v n+1/2 At (38) 

V n +1 = v n+l/2 + f(x n+ i)Z\i/2 (39) 
or a DKD version of the method 

x n+1/2 =x„ + v„^t/2 (40) 

v„+i =v„+f(x„ +1/2 )zii (41) 



x 



n+1 = x, i+ i /2 + v n+i At/2. (42) 

This method is accurate to second order, and, as will be shown in the next paragraph, 
also has other advantages. For more details see ISpringell l|2005l ). 



It is also clear that, depending on the application, a lower order scheme applied with 
more, and thus smaller, time steps can be more efficient than a higher order scheme, 
which enables the use of larger time steps. In the upper rows of Fig. |4j we show the 
numerical integration of a Kepler problem (i.e. two point-like masses with large mass 
difference which orbit around each other like a planet-sun system) of high eccentricity 
e = 0.9, using second-order accurate leap-frog and Runga-Kutta schemes with fixed 
time step. There is no long-term drift in the orbital energy for the leap-frog result (left 
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panel); only a small residual precession of the elliptical orbit is observed. On the other 
hand, the second-order Runge-Kutta integrator, which has formally the same error per 
step, fails catastrophically for an equally large time step (middle panel). After only 50 
orbits, the binding energy has increased by ~ 30 %. If we instead employ a fourth-order 
Runge-Kutta scheme using the same time step (right panel), the integration is only 
marginally more stable, now giving a decline of the binding energy of ~ 40 % over 
200 orbits. Note however that such a higher order integration scheme requires several 
force evaluations per time step, making it computationally much more expensive for a 
single step than the leap-frog, which requires only one force evaluation per step. The 
underlying mathematical reason for the remarkable stability of the l eap-frog integra tor 
lies in its symplectic properties. For a more detailed discussion, see ISpringell (2005). 

In cosmological simulations, we are confronted with a large dynamic range in 
timescales. In high-density regions, like at the centres of galaxies, the required time 
steps are orders of magnitude smaller than in the low-density regions of the inter- 
galactic medium, where a large fraction of the mass resides. Hence, evolving all the 
particles with the smallest required time step implies a substantial waste of computa- 
tional resources. An integration scheme with individual time steps tries to cope with 
this situation more efficiently. The principal idea is to compute forces only for a certain 
group of particles in a given kick operation (K), with the other particles being evolved 
on larger time steps being usually just drifted (D) and 'kicked' more rarely. 

The KDK scheme is hence clearly superior once one allows for individual time 
steps, as shown in the lower ro w of Fig. [4| It is also possible to try to recover the 
time reversibility more precisely. iHut et aL (|l995h discuss an implicit time step cri- 
terion th at depends both on the beginning and on the end of the time step, and, 
similarly, iQuinn et al.l ( 19971 ) discuss a binary hierarchy of trial steps that serves a 
similar purpose. However, these schemes are computationally impractical for large col- 
lisionless systems. Fortunately, however, in this case, the danger of building up large 
errors by systematic accumulation over many periodic orbits is much smaller, because 
the gravitational potential is highly time-dependent and the particles tend to make 
comparatively few orbits over a Hubble time. 



2.6 Initial conditions 

Having robust and well justified initial conditions is one of the key points of any 
numerical effort. For cosmological purposes, observations of the large-scale distribution 
of galaxies and of the CMB agree to good precision with the theoretical expectation 
that the growth of structures starts from a Gaussian random field of initial density 
fluctuations; this field is thus completely described by the power spectrum P(|k|) whose 
shape is theoretically well motivated and depends on the cosmological parameters and 
on the nature of Dark Matter. 

To generate the initial conditions, one has to generate a set of complex numbers 
with a randomly distributed phase <f) and with amplitude norm ally distributed with a 
variance given by the desired spectrum (e.g. lBardeen et al . 1986). This can be obtained 
by drawing two random numbers <f> in ]0, 1] and A in ]0, 1] for every point in fe-space 



4 = v/-2P(|k|)hi(4)e l2 



(43) 
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Fig. 5 Shown is a slice to the particle distribution with the imposed displacement, taken from 
the same cosmological initial conditions, once based on an originally regular grid (left panel) 
and once based on an originally glass like particle distribution (right panel). 



To obtain the perturbation field generated from this distribution, one needs to generate 
the potential <P(q) on a grid q in real space via a Fourier transform, e.g. 



*(q)=£^ 



(44) 



The subsequent application of the Zel'dovich approximation ( Zerdovichlll970l ) enables 
one to find the initial positions 



and velocities 



x = q-Z? + (z)<£(q) 



v = D + {z)V<P(q) 



(45) 



(46) 



of the particles, where D + (z) and D + (z) indicate the cosmological linear growth factor 
and its derivative at t he ini tial redshift z. A more detailed description can be found in 
e.g. lEfstathiou et al.l (|l985l ). 

There are two further complications which should be mentioned. The first is that 
one can try to reduce the discreteness effect that is induced on the density power 
spectrum by the regularity of the underlying grid of the particle positions q that one 
has at the start. This can be done by constructing an amorphous, fully relaxed particle 
distribution to be used, instead of a regular grid. Such a particle distribution can be 
constructed by applying negative gravity to a system and evolving it for a long time, 
including a d amping of the velocities, until it reaches a relaxed state, as suggested by 
IWhiteli|l99ri ). Fig. © gives a visual impression on the resulting particle distributions. 

A second complication is that, even for studying individual objects like galaxy 
clusters, large-scale tidal forces can be important. A common approach used to deal 
with this problem is the so-called "zoom" technique: a high resolution region is self- 
consistently embed ded in a larger scale cosmological volume at low resolution (see e.g. 
iTormen et ail 19971 ). This approach usually allows an increase of the dynamical range of 
one to two orders of magnitude while keeping the full cosmological context. For galaxy 
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Fig. 6 Mean inner density contrast as a function of the enclosed number of particles in 4 series 
of simulations varying the number of particles in the high-resolution box, from 32 3 to 256 3 . 
Each symbol corresponds to a fixed fraction of the virial radius, as shown by the labels on the 
right. The number of particles needed to obtain robust results increases with density contrast, 
roughly as prescribed by the requirement that the collisional relaxation timescale should remain 
longer than the age of the Universe. According to this, robust numerical estimates of the mass 
profile of a halo ar e only possible to the right of the curve labelled t rc lax ~ 0.64o- Taken from 
iPower et all l|2003l 1. 



simulations it is even possible to apply this technique on sever al levels of refinem ents 
to further improve the dynamical range of the simulation (e.g. IStoehr et al.ll2003h . A 
frequently used, publicly availa ble package to create initial conditions is the COSMICS 
package bv lBertschirigerl (jl995f l. 



2.7 Resolution 

There has been a long standing discussion in the literature to understand what is 
the optimal setup for cosmological simulations, and how many particles are needed 
to resolve certain regions of interest. Note that the number of particles needed for 
convergence also depends on what quantity one is interested in. For example, mass 
functions, which count identified haloes, usually give converging results at very small 
particle numbers per halo (~ 30 — 50), whereas structural properties, like a central 
density or the virial radius, converge only at significantly higher particle numbers 
(~ 1000). As we will see in a later chapter, if one wants to infer hydrodynamical 
properties like baryon fraction or X-ray luminosity, values converge only for haloes 
represented by even more partic les (~ 10 000). 

Recently. Fower et al.l (|200&t ) performed a comprehensive series of convergence tests 
designed to study the effect of numerical parameters on the structure of simulated 
CDM haloes. These tests explore the influence of the gravitational softening, the time 
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Fig. 7 Moore's empirical law shows that the computing power typically doubles every 18 
months. This figure shows the size of N-body simulations as a function of their running date. 
Clearly, specially recently, the improvement in the algorithms allowed the simulation to grow 
faster than the improvement of the underlying CPU power. Kindly provided by Volker Springel. 



stepping algorithm, the starting redshift, the accuracy of force computations, and the 
number of particles in the spherically-averaged mass profile of a gal axy-sized hal o in th e 
CDM cosmogony with a non-null cosmological constant fylCDM). IPower et al.l (120051 ) . 
and the references therein, suggest empirical rules that optimise the ch oice of these 
parameters. When these choices arc dictated by computational limitations. Power et al. 
(2003) offer simple prescriptions to assess the effective convergence of the mass profile 
of a simulated halo. One of their main results is summarised in Fig. [6] which shows 
the convergence of a series of simulations with different mass resolution on different 
parts of the density profile of a collapsed object. This figure clearly demonstrates that 
the number of particles within a certain radius needed to obtain converging results 
depends on the enclosed density. 

In general, both the size and the dynamical range or resolution of the simulations 
have been increasing very rapidly over the last decades. Fig. [7] shows a historical com- 
pilation of large N-body simulations: their size growth, thanks to improvements in the 
algorithms, is faster than the underlying growth of the available CPU power. 
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Fig. 8 A recent comparison of the predicted number of halos as a function of density for ten 
different cosmological codes. Left panel: halos with 10 — 40 particles, right panel: halos with 
41 — 2500 particles. The lower panels show the residuals with respect to GADGET-2. Both 
panels show the deficit of small halos in ENZO and FLASH over most of the density region - 
only at very high densities do the results catch up. The behaviour of the TPM simulation is 
interesting: not only does this simulation have a deficit of small halos but the deficit is very 
significant in medium density regions, in fact falling below the two Adaptive Mesh Refinement 
codes. The slight excess of small halos shown in the TreePM run vanishes completely if the 
halo cut is raised to 20 particles per halo an d the TreePM re s ults a re in that case in excellent 
agreement with GADGET-2. Adapted from iHeitmann et all 1120071) . 



2.8 Code comparison for pure gravity 

In the last thirty years cosmology has turned from a science of order-of-magnitude 
estimates to a science with accuracies of 10 % or less in its measurements and theo- 
retical predictions. Crucial observations along the way were the measurement of the 
cosmic microwave background radiation, and large galaxy surveys. In the future such 
observations will yield even higher accuracy (1 %) over a much wider dynamical range. 
Such measurements will provide insight into several topics, e.g. the nature of dark en- 
ergy (expressed by the equation of state w = p/p with p being the pressure and p the 
density). In order to make optimal use of the observations, theoretical calculations of 
at least the same level of accuracy are required. As physics in the highly non-linear 
regime, combined with complicated gas physics and astrophysical feedback processes 
are involved, this represents a real challenge. 

Different numerical methods have therefore to be checked and compared contin- 
uously. The most recent comparison of ten commonly-used codes from the literature 
has been performed in a n extensive comparison program. The ten codes used for the 
comparison performed by IHeitmann et all ( 20071 ) cover a variety of methods and appli- 



cation arenas. The simulation methods employed include parallel particle-in-cell (PIC) 
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techniques (the PM codes MC 2 and PMM, the Particle-Mesh / Adaptive Mesh Re- 
finement (AMR) codes ENZO and FLASH), a hybrid of PIC and direct N-body (the 
AP 3 M code Hydra), tree algorithms (the treecodes PKDGRAV and HOT), and hybrid 
tree-PM algorithms (GADGET-2, TPM, and TreePM). 

The results from the code comparisons are satisfactory and not unexpected, but 
also show that much more work is needed in order to attain the required accuracy 
for upcoming surveys. The halo mass function is a very stable statistic, the agreement 
over wide ranges of mass being better than 5 %. Additionally, the low mass cutoff for 
individual codes can be reliably predicted by a simple criterion. 

The internal structure of halos in the outer regions of ~ -R200 a l so appears to be 
very similar between different simulation codes. Larger differences between the codes in 
the inner region of the halos occur if the halo is not in a relaxed state: in this case, time 
stepping issues might also play an important role (e.g. particle orbit phase errors, global 
time mismatches). For halos with a clear single centre, the agreement is very good and 
predictions for the fall-off of the profiles from resolution criteria hold as expected. The 
investigation of the halo counts as a function of density revealed an interesting problem 
with the TPM code, the simulation suffering from a large deficit in medium density 
regimes. The AMR codes showed a large deficit of small halos over almost the entire 
density regime, as the base grid of the AMR simulation sets a resolution limit that is 
too low for the halos, as can be seen in Figure ((8}. 

The power spectrum measurements revealed definitively more scatter among the 
different codes than expected. The agreement in the nonlinear regime is at the 5 — fO % 
level, even on moderate spatial scales around k — Wh Mpc -1 . This disagreement on 
small scales is connected to diff erences of the codes in the inner regions of the halos. 
For more detailed discussion see Hcitm ann et al. and references therein. 

In a detailed comparison of ENZO and GADGET, lO'Shea et all (|2005h already 
pointed out that to reach reasonable good agreement, relatively conservative criteria for 
the adaptive grid refinement are needed. Furthermore, choosing a grid resolution twice 
as high as the mean inter-particle distance of the dark matter particles is recommended, 
to improve the small scale accuracy of the calculation of the gravitational forces. 



3 Hydro methods 

The baryonic content of the Universe can typically be described as an ideal fluid. 
Therefore, to follow the evolution of the fluid, one usually has to solve the set of 
hydrodynamic equations 

dv VP 

— = V<P, 47 

dt p 

-£ + pVv = (48) 

and 

^=_£ V .v-^), (49) 
dt p p v ; 

which are the Euler equation, continuity equation and the first law of thermodynamics, 
respectively. They are closed by an equation of state, relating the pressure P to the 
internal energy (per unit mass) u. Assuming an ideal, monatomic gas, this will be 



P = (7 - l)pu 



(50) 



19 



f „,„(*) = 



PCM 
PLM 
PPM 



u n-2 



u n-1.5 u n-l 



u n+O.S 




u n+1.5 "n+2 



u „ + l u n + l.S 



Fig. 9 Reconstruction of the principal variables on the grid using different methods, 

like piecewise constant (PCM), piecewise linear (PLM) or piecewise parabolic (PPM). The 
reconstruction scheme then allows one to calculate cell averages (u n ) as well as the left and 
right-hand sided values on the cell boundaries (mJ^q 51^-1-0 5)' 



with 7 = 5/3. In the next sections, we will discuss how to solve this set of equations, 
neglecting radiative losses described by the cooling function A(u, p); in Sect. 4.1 we 
will give examples of how radiative losses or additional sources of heat are included 
in cosmological codes. We can also assume that the V# term will be solved using the 
methods described in the previous section. 

As a result of the high nonlinearity of gravitational clustering in the Universe, 
there are two significant features emerging in cosmological hydrodynamic flows; these 
features pose more challenges than the typical hydrodynamic simulation without self- 
gravity. One significant feature is the extremely supersonic motion around the density 
peaks developed by gravitational instability, which leads to strong shock discontinuities 
within complex smooth structures. Another feature is the appearance of an enormous 
dynamic range in space and time, as well as in the related gas quantities. For instance, 
the hierarchical structures in the galaxy distribution span a wide range of length scales, 
from the few kiloparsecs resolved in an individual galaxy to the several tens of mega- 
parsecs characterising the largest coherent scale in the Universe. 

A variety of numerical schemes for solving the coupled system of collisional baryonic 
matter and collisionless dark matter have been developed in the past decades. They fall 
into two categories: particle methods, which discretise mass, and grid-based methods, 
which discretise space. We will briefly describe both methods in the next two sections. 
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and 



3.1 Eulerian (grid) 

The set of hydro-dynamical equations for an expanding Universe reads 
<9v 1 , a 1 __ 1_. 

— - + -(v ■ V)v + -v = VP V<5, 51) 

at a a ap a 

g + ^+V(pv) = (52) 

-|(H + ~v ■ V(pu) = -(pu + P) (iv ■ v + 3^) (53) 

respectively, where the right term in the last equation reflects the expansion in addition 
to the usual PdV work. 

The grid-based methods solve these equations based on structured or unstructured 
grids, representing the fluid. One distinguishes primitive variables, which determine 
the thermodynamic properties, (e.g p, v or P) and conservative variables which define 
the conservation laws, (e.g. p, pv or pu). Early attempts were made using a central 
difference scheme, where fluid is only represented by the centred cell values (e.g. central 
variables, u n in Fig. [9] and derivatives are obt ained by th e finite-difference representa- 
tion, similar to Eg. I15l and ll6l see for example ICenl l|l992h . Such methods will however 



break down in regimes where discontinuities appear. These methods therefore use ar- 
tificial viscosity to handle shocks (similar to the smoothed particle hydrodynamics 
method described in the next section). Also, by construction, they are only first-order 
accurate. 

More modern approaches use reconstruction schemes, which, depending on their 
order, take several neighbouring cells into account to reconstruct the field of any hy- 
drodynamical variable. Fig. [9] illustrates three different reconstruction schemes, with 
increasing order of ac curacy, as piecewise constant method (PCM), piecewise lin- 
ear method (PLM, e.g. Colella fe Glad[l985l ) and piecewise parabolic method (PPM, 



IColella fc Woodward! il984T T The shape of the reconstruction f n , u (x) is then used to 



calculate the total integral of a quantity over the grid cell, divided by the volume 
of each cell (e.g. cell average, u n ), rather than pointwise approximations at the grid 
centres (e.g. central variables, u n ). 

2„ + 0.5 

/n, u (a;)da; (54) 



They are also used to calculate the left and right-hand sided values at the cell 
boundaries (e.g. u l n ± 5 ,it^ ±0 5 ), which are used later as initial conditions to solve the 
Riemann problem. To avoid oscillations (e.g. the development of new extrema), addi- 
tional constraints are included in the reconstruction. For example, in the PLM recon- 
struction this is ensured by using so-called slope limiters which estimate the maximum 
slope allowed for the reconstruction. One way is to demand that the total variation 
among the interfaces does not incre ase with time. Such so-called total variation dimin- 
ishing schemes ( TVD, iHartenll 19831 ') . nowadays provide various different slope limiters 
suggested by different authors. In our example illustrated by Fig. [5] the so called min- 
mod slope limiter 

Aui = minmod (Q(u i+1 - Uj), (u i+1 - it i _ 1 )/2, - Wi-i)) , (55) 
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where Azij is the limiter slope within the cell i and 9 — [1,2], would try to fix the slope 
f'n—X u( x n—l) an d fn,u( x n), such as to avoid that u l n _ 5 becomes larger than u^_ 5- 
The so called Aldaba-type limiter 

where e is a small positive number to avoid problems in homogeneous regions, would 
try to avoid that u l n _ 5 is getting larger than u n and that u! n _ 5 is getting smaller 
than U n _i, e.g. that a monotonic profile in Ui is preserved. 

In the PPM (or even higher order) reconstruction this enters as an additional con- 
dition when finding the best-fitting polynomial function. The additional cells which 
are involved in the reconstruction are often called the stencil. Modern, high order 
schemes usually have stenci ls based on at leas t 5 grid points and implement essen- 
tially non-oscillatory (ENO; lHarten et al.| [l987) or monoticity preserving (MP) meth- 
ods for reconstruction, which maintain high-order accuracy. For every reconstruction, 
a smoothness indicator S™ can be constructed, which is defined as the integral over 
the sum of the squared derivatives of the reconstruction over the stencil chosen, e.g. 

2 x n + m 

S™=Y1 I i^f- 1 (difnA*)) 2 (57) 

' =1 x„ m 

In the ENO schemes, a set of candidate polynomials with order 2m + 1 for a set of 
stencils based on different numbers of grid cells m are used to define several different 
reconstruction functions /,'"„. Then, the reconstruction with the lowest smoothness 
indicator S 1 ™ is chosen. In this way the order of reconstruction will be reduced around 
discontinuities, and oscillating behaviour will be suppressed. 

To improve on the ENO schemes in robustness and accuracy one can, instead of 
selecting the reconstruction with the best smoothness indicator S™, construct the final 
reconstruction by building the weighted reconstruction 

fu(x) = 2J Wmfn,u{x), (58) 
m 

where the weights w m are a proper functio n of the smoothness indicators S™. This 
procedure is not unique, [jiang fc Shul (1996) proposed defining 



with 



w m = J^ m (59) 



<*i = %t^, (60) 



where Q, e and /3 are free parameters, which for example can be taken from lLevv et al.l 
(1999). This are the so-called weighted essentially non-oscillatory (WENO) schemes. 



These schemes can simultaneously provide a high-order resolution for the smooth part 
of the solution and a sharp, monotonic shock or contact d iscontinuity transition. For 
a review on ENO and WENO schemes, see e.g. Ishul l| 19981 ). 

After the left and right-hand values at the cell boundaries (e.g. interfaces) are re- 
constructed, the resulting Riemann problem is solved, e.g. the evolution of two constant 
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Fig. 10 The Ricmann problem: The upper panel shows the initial state, the lower panel shows 
the evolved problem for the case of no relative motion between the two sides (ui = u$ = 0). 
The solid lines mark the pressure P, the dashed dotted lines the density p and the dotted line 
the velocity v. Kindly provided by Ewald Miiller. 



states separated by a discontinuity. This can be done either analytically or approxi- 
mately, using left and right-handed values at the interfaces as a jump condition. 

With the solution one obtains, the fluxes across these boundaries for the time step 
can be calculated and the cell averages u n can be updated accordingly. In multiple di- 
mensions, all these steps are performed for each coordinate direction separately, taking 
only the interface values along the individual axes into account. There are attempts 
to extend the reconstruction schemes, to directly reconstruct the principal axis of the 
Riemann problem in multiple dimensions, so that then it has to be solved only once 
for each cell. However the complexity of reconstructing the surface of shocks in three 
dimensions has so far seen to be untraceable. 

How to solve the general Riemann problem, e.g. the evolut ion of a discontinuity 
initia lly separating two states, can be found in text books (e.g. ICourant fc Friedrichj 
1 19481 ). Here we give only the evolution of a shock tube as an example. This corresponds 
to a system where both sides are initially at rest. Fig. [TO] shows the initial and the 
evolved system. The latter can be divided into 5 regions. The values for regions 1 and 
5 are identical to the initial configuration. Region 2 is a rarefaction wave which is 
determined by the states in region 1 and 3. Therefore we are left with 6 variables to 
determine, namely P3, P3, 113 and p4, P4, V4, where we have already eliminated the 
internal energy Uj in all regions, as it can be calculated from the equation of state. 
As there is no mass flux through the contact discontinuity, and as the pressure is 
continuous across the contact discontinuity, we can eliminate two of the six variables 
by setting v$ = 774 = v c and P3 = P4 = P c . The general Rankine-Hugoniot conditions, 
describing the jump conditions at a discontinuity, read 



PlVi = p r v r 
Pm + Pi = PrV r + Pr 



(61) 

(62) 
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VI (pi{vf/2 + U{) + Pi) = Vr ( P r(vr/2 + U r ) + Pr) , (63) 

where we have assumed a coordinate system which moves with the shock velocity v s . 
Assuming that the system is at rest in the beginning, e.g. ^1=^5= 0, the first Rankine- 
Hugoniot condition for the shock between region 4 and 5 moving with a velocity v s 
(note the implied change of the coordinate system) is in our case 

m = P5 v s — pi(v s — v c ) (64) 
and therefore the shock velocity becomes 

P4Vc , _s 

v s = • (65) 

Pi - P5 

The second Rankine-Hugoniot condition is 

mvc = pi(v s — v c )vc = Pc — P5, (66) 
which, combined with the first, can be written as 

p 4 f^^_ _ v \ Vc = Pc - p 5 , ( 6 7) 
\P4-P5 J 

which, slightly simplified, leads to a first condition 

(P 5 -Pc)(---) =-v 2 c . (68) 

\P5 Pi J 



The third Rankine-Hugoniot condition is 

2 



m(6 4 + ^ - £ 5 ) = Pcvc, (69) 
which, by eliminating m, can be written as 

Pc + P5 Vc f7n , 

e ^ e5= p7^p- 5 Y- (70) 

Using the first condition (Eq. I68[) and assuming an ideal gas for the equation of state, 
one gets 

1 f Pc P 5 \ Pc + P 5 



7 — 1 \_ P4 P5 J 2p4P5 

which leads to the second condition 



(71) 



Pc " P5 =7 P4~P5 (72) 
Pc+P 5 'P4 + P5 ^ ' 

The third condition comes from the fact that the entropy (oc ln(P/p 7 ) stays constant 
in the rarefaction wave, and therefore one can write it as 



Pi _ ( Pi X 



Pc \ P3/ 

The fourth condition comes from the fact that the Riemann Invariant 



(73) 



r + I (74) 
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is a constant, which means that 



v c + I — dp = vi+ —dp, 
P3 J PI 



(75) 



where c = \fyP/p denotes the sound velocity, with which the integral can be written 

(76) 



Therefore, the fourth condition can be written as 



v c + 



1 - 1 V P3 



7-1 



iPl 
Pi 



(77) 



Combining all 4 conditions (eauations l68ll72ll73l and l77[l and defining the initial density 
ratio A = P1/P5 one gets the non linear, algebraic equation 



pi i (i-py 



2 7 



P5 A 7 (l + P)-1 + P (7-I) 2 



p\ (7"l)/(27) 

aJ 



(78) 



for the pressure ratio P = P c /P§- Once P c is known from solving this equation, the 
remaining unknowns can be inferred step by step from the four conditions. 

There are various approximate methods to solve the Riemann problem, includ- 
ing the so-called ROE method (e.g. IPowell et alj|l999h . HLL/HLLE method (e.g. see 



lHarten et allll983i : lEinfeid3ll988i ; lEinfeldt et al.lll99ll ) and HLLC (e.g. see 111 120051 ) ■ A 

description of all these methods is outside the scope of this rev iew, so we redirect the 
reader to the references given or textbooks like iLeVeauel (|2002T ). 

At the end of each time step, one has to compute the updated central values u n from 
the updated cell average values u n . Normally, this would imply inverting equation (|54l) . 
which is not trivial in the general case. Therefore, usually an additional constraint is 
placed on the reconstruction method, namely that the reconstruction fulfills u n — u n . 
In this case the last step is trivial. 

In general, the grid-based methods suffer from limited spatial resolution, but they 
work extremely well in both low- and high-density regions, as well as in shocks. In 
cosmological simulations, accretion flows with large Mach numbers (e.g. M > 100) are 
very common. Here, following the total energy in the hydrodynamical equations, one 
can get inaccurate thermal energy, leading to negative pressure, due to discretisation 
errors when the k ineti c energy dominates the total energy. In such cases, as suggested by 
iRvu et all 1 1993h and lBrvan et alj (|l995h . the numerical schemes usually switch from 
formulations solving the total energy to formulations based on solving the internal 
energy in these hypersonic flow regions. 

In the cosmo logic al setting , ther e are the TVD-based codes including tho se of 
Rvu et al.1 (|l993h and iLi et aTT (|2006 ) (CosmoMHD) the moving-mesh s cheme (|Pe 



19981 ) and the PLM-based c ode ART (iKravtsoy et"allll997l ; lKravtsovll2002h. Th e PPM- 
based codes include those o f lStone fc Normanldl992l ) (Zeu s) lBrvan et al.l(|l995h (ENZO), 



iRicker et all |2000h ( COSMOS) andlFryxell et all (|2000h (FLASH). There is also the 
WEJVO-based code bv lFeng et ail (|2004l ). 
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3.2 Langrangian (SPH) 



The particle methods include variant s of smoothed particle hydr o dynamics (SPH; 
iGineold fc Monaghanll977 :ILucvIi977I) such as those oflEvrardl Jl98Sl).lHernauist fc Katzl 
( 19891 ), iNavarro fc White! Jl993h. ICouchman et"afl dl995h ("Hy dra), Istejnmetj |l996g~ 



Owen et al 



( 19981), and lSpringel et all (|2001ah : ISpringei l|200a i (GAD- 



(GRAPESPH), 

GET). The SPH method solves the Lagrangian form of the Euler equations and can 
achieve good spatial resolutions in high-density regions, but it works poorly in low- 
density regions. It also suffers from degrad ed resolution in sho cked regions due to the 
introduction of a sizable artificial viscosity. lAgertz et all (|2007h argued that whilst Eu- 
lerian grid-based methods are able to resolve and treat dynamical instabilities, such as 
Kelvin-Helmholtz or Rayleigh-Taylor, these processes are poorly resolved by existing 
SPH techniques. The reason for this is that SPH, at least in its standard implementa- 
tion, introduces spurious pressure forces on particles in regions where there are steep 
density gradients, in particular near contact discontinuities. This results in a boundary 
gap of the size of an SPH smoothing kernel radius, over which interactions are severely 
damped. Nevertheless, in the cosmological context, the adaptive nature of the SPH 
method compensates for such shortcomings, thus making SPH the most commonly 
used method in numerical hydrodynamical cosmology. 



3.2.1 Basics of SPH 

The basic idea of SPH is to discretise the fluid by mass elements (e.g. particles), rather 
than by volume elements as in the Eulerian methods. Therefore it is immediately 
clear that the mean inter-particle distance in collapsed objects will be smaller than in 
underdense regions; the scheme will thus be adaptive in s patial resolution b y keeping 
the mass resolution fixed. For a comprehensive review see iMonaghanl l|l992l ). To build 



continuous fluid quantities, one starts with a general definition of a kernel smoothing 
method 

(A(x)) = / W(x - x', h)A(x!)dx, (79) 



which requires that the kernel is normalised (i.e. J VK(x, h)dx — 1) and collapses to a 
delta function if the smoothing length h approaches zero, namely W(x.,h) — > <5(x) for 
h -> 0. 

One can write down the continuous fluid quantities (e.g. (A(x)}) based on the 
discretised values Aj represented by the set of the individual particles rrij at the position 

(A,) = (A( Xi )) = V — AjW(xj - xj, h) , (80) 

where we assume that the kernel depends only on the distance modulus (i.e. W(|x — 
x'\,h)) and we replace the volume element of the integration, dx = d 3 x, with the 
ratio of the mass and density mj/pj of the particles. Although this equation holds 
for any position x in space, here we are only interested in the fluid representation at 
the original particle positions Xj, which are the only locations where we will need the 
fluid representation later on. It is important to note that for kernels with compact 
support (i.e. W(x,h) = for |x| > h) the summation does not have to be done over 
all the particles, but only over the particles within the sphere of radius h, namely the 
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neighbours around the particle i under consideration. Traditionally, the most frequently 
used kernel is the i?2-sphne, which can be written as 

a [ i-Htf + Hrf o<£<o.5 

W{x,h) = ^\2(l-%f 0.5<f<l> (81) 

where v is the dimensionality (e.g. 1, 2 or 3) and a is the normalisation 

x 4 » = i 

§ u = 2 (82) 
I - = 3. 

Sometimes, spline kernels of higher order are used for very special applications; however 
the £?2 spline kernel turns out to be the optimal choice in most cases. 

When one identifies Aj with the density p$, pi cancels out on the right hand side 
of Eq. 1801 and we are left with the density estimate 

{p i )=J2m j W(x i -x i ,h), (83) 



which we can interpret as the density of the fluid element represented by the particle 
i. 

Now even derivatives can be calculated as 

V{Ai)=y2—A j V i W(x l —x i ,h), (84) 
i Pj 

where V» denotes the derivative with respect to x,. A pairwise symmetric formulation 
of derivatives in SPH can be obtained by making use of the identity 

(pV)-A = V(p-A)-p-(VA), (85) 

which allows one to re-write a derivative as 

V (A) = - V m 3 (A 3 - Ai)ViW(xi - xj, h). (86) 

3 

Another way of symmetrising the derivative is to use the identity 

^=vfi) + 4v ft (87) 
which then leads to the following form of the derivative: 

V{A i )=p i Y^m j (^yr + ^J ViW(xi-xj,A). (88) 
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3.2.2 The fluid equations 



By making use of these identities, the Euler equation can be written as 

^ = -E^R + 5 + ^'] ViW(Xi - x j- (89) 



By combining the above identities and averaging the result, the term — (P/p)V ■ v from 
the first law of thermodynamics can similarly be written as 

-v^ViWCxi-xj.Ti). (90) 

Here we have added a term 77^ which is the so-called artificial viscosity. This term 
is usually needed to capture shocks and its construction is si milar to other hydro- 
dynam ical s chemes. Usuall y, one adopts the form proposed by iMonaghan &i Gingoldl 
( 19831 ) and iBalsaral (|l995l ), which includes a bulk viscosity and a von Neumann- 




Richtmeyer viscosity term, supplemented by a term controlling a ngular momentum 
transport in the presence of shear flows at low particle numbers |Steinme tz||l996bl) ■ 



Mode rn schemes implement a form of the artificial viscosity as proposed by MonagharJ 



( 19971 ) aased on an analogy with Riemann solutions of compressible gas dynamics. To 



reduce this artificial viscosity, at least in tho se parts of the flows where there are no 
shocks, one can follow the idea proposed bv iMorris fc MonagharJ 1 1997t ): every par- 



ticle carries its own artificial viscosity, which eventually decays outside the regions 
which undergo shocks. A detail ed study of t he im plications on the ICM of such an 
implementation can be found in iDolag et al.l (|2005h . 



The continuity equation does not have to be evolved explicitly, as it is automatically 
fulfilled in Lagrangian methods. As shown earlier, density is no longer a variable but can 
be, at any point, calculated from the particle positions. Obviously, mass conservation is 
guaranteed, unlike volume conservation: in other words, the sum of the volume elements 
associated with all of the particles might vary with time, especially when strong density 
gradients are present. 



3.2.3 Variable smoothing length 



Usually, the smoothing length h will be allowed to vary for each individual particle i 
and is determined by finding the radius hi of a sphere which contains n neighbours. 
Typically, different numbers n of neighbours are chosen by different authors, ranging 
from 32 to 80. In p rinciple, dependi ng on the kernel, there is an optimal choice of 
neighbours (e.g. see ISilverman] Il98fj| or similar books) . However, one has to find a 
compromise between a large number of neighbours, leading to larger systematics but 
lower noise in the density estimates (especially in regions with large density gradients) 
and a small number of neighbours, leading to larger sample variances for the density 
estimation. In general, once every particle has its own smoothing length, a symmetric 
kernel W(x.i — Xj, hi, hj) — Wij has to be constructed to keep the conservative form of 
the formulations of the hydrodynamical equations. There are two main variants used 
in the literature: one is the kernel average Wij = (W^Xj — x^-, hi) + W(x.i — Xj, hj))/2, 
the other is an average of the smoothing length Wij = W(~x.i — Xj, (hi + hj)/2). The 
former is the most commonly used approach. 
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Note that in all of the derivatives discussed above, it is assumed that h does not 
depend on the position Xj . Thus, by allowing the smoothing length hi to be variable for 
each particle, one formally neglects the correction term dW/dh, which would appear 
in all the derivatives. In general, this correction term cannot be computed trivially and 
therefore many implementations do not take it into account. It is well known that such 
formulations are poor at conserving numerically both internal energy and entropy at 
the same time, independently of the use of internal energy or entropy in the formulation 
of the first law of thermodynamics, see iHernauistJ (|l993h . In the next subsection, we 
present a way of deriving the equations which include these correction terms dW/dh; 
this equation set represents a formulation which conserves numerically both entropy 
and internal energy. 



3.2.4 The entropy conservation formalism 

To derive a better formulation of the SPH method, ISpringel fc Hernquistl (|2002h started 
from the entropic function A = Pip 1 , which will be conserved in adiabatic flows. The 
internal energy per unit mass can be inferred from this entropic function as 

X% TPV (91) 



7- 

at any time, if needed. Entropy will be generated by shocks, which are captured by the 
artificial viscosity 77^ and therefore the entropic function will evolve as 

IF = E m i n v ( V J - v v *t% ( 92 ) 

The Euler equation can be derived starting by defining the Lagrangian of the fluid as 

£(q, q) = 2 mii * — ~~i X! m^Pi' 1 (93) 

i i 

which represents the entire fluid and has the coordinates q = (xi, x^r, hi, hpj). 
The next important step is to define constraints, which allow an unambiguous associ- 
ation of hi for a chosen number of neighbours n. This can be done by requiring that 
the kernel volume contains a constant mass for the estimated density, 

47T -J 

&(q) = yfciPi - nm * = °- (94) 



dt dq, den j den 



The equation of motion can be obtained as the solution of 

which - as demonstrated bv lSpringel fc Hernquistl (|2002h - can be written as 
dv, 



3 



— - - / nvj . jj— V l W(x i -Kj,hj) + fi— ViW(x; - X j , hi) + IlijV \Wij 



(96) 
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Fig. 11 One-dimensional velocity dispersion profile (left panel) and gas temper ature profile 
(righ t panel) of the cluster at z = of the Santa Barbara Comparison Project l lFrenk et al.l 
1999). The solid line is the profile averaged over the 12 simulations. The symbols correspond 
to individual simulations. The crosses in the left panel correspond to a dark- matter only 
simula tion. The top panels show the residual from the mean profile. Taken from lFrenk et all 
41999T ). 



where we already have included the additional term due to the artificial viscosity 77^ , 
which is needed to capture shocks. The coefficients /.; incorporate fully the variable 
smoothing length correction term and are defined as 

Note that in addition to the correction terms, which can be easily calculated together 
with the density estimation, this formalism also avoids all the ambiguities we saw in 
the derivations of the equations in the previous section. This formalism defines how the 
kernel averages (symmetrisation) have to be taken, and also fixes hi to unambiguous 
values. For a detailed deriv ation of this formalism and its conserving capabilities see 
ISpringel fc Hernquistl |2002l ). 



3.3 Code Comparison 

The Eulerian and Lagrangian approaches described in the previous sections should 
provide the same r esults when applie d to the same problem, like the interaction of 
multi-phase fluids ( Agertz et a.1.1 [2007^ . To verify that the code correctly solves the 
hydrodynamical set of equations, each code is usually tested against problems whose 
solution is known analytically. In practice, these are shock tubes or spherical collapse 
problems. In cosmology, a relevant test is to compare the results provided by the 
codes when they simulate the formation of cosmic structure, when finding an analytic 
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solution is impractical; for example lO'Shea et al.l 1 2005h compare the thermodynamical 
properties of the intergalactic medium predicted by the GADGET (SPH-based) and 
ENZO (grid-based) codes. Anoth er example of a co mparison between grid-based and 
SPH-based codes can be found in lKang et al. I 1 19941 ). 

A detailed comparison of hydrodynamical codes which simulate the formation and 
evolution of a clus t er was provided by the Santa Barbara Cluster Comparison Project 
( Frenk et al.lll999l ). lFrenk et alJ (1999) comprised 12 different groups, each using a code 
either based on the SPH technique (7 groups) or on the grid technique (5 groups). Each 
simulation started with identical initial conditions of an individual massive cluster in a 
flat CDM model with zero cosmological constant. Each group was free to decide resolu- 
tion, boundary conditions and the other free parameters of their code. The simulations 
were performed ignoring radiative losses and the simulated clusters were compared at 
z = 0.5 and z = 0. 

The resulting dark matter properties were similar: it was found a 20 % scatter 
around the mean density and velocity dispersion profiles (left panel of Fig. Ill[) . A 
similar agreement was also obtained for many of the gas properties, like the temperature 
profile (right panel of Fig. Hip or the ratio of the specific dark matter kinetic energy 
and the gas thermal energy. 

Somewhat larger differences are present for the inner part of the temperature or 
entropy profiles and more recent implementations have not yet cured this problem. The 
largest discrepancy was in the total X-ray luminosity. This quantity is proportional 
to the square of the gas density, and resolving the cluster central region within the 
core radius is crucial: the simulations resolving this region had a spread of 2.6 in the 
total X-r ay luminosit y , com pared to a spread of 10 when all the simulations were 
included. iFrenk et al.1 (11999) also concluded that a large fraction of the discrepancy, 
when excluding the X-ray luminosity result, was due to differences in the internal 
timing of the simulations: these differences produce artificial time shifts between the 
outputs of the various simulations even if the outputs are formally at the same cosmic 
time. This reflects mainly the underlying dark matter treatment, including chosen force 
accuracy, different integration schemes and choice of time steps used, as described 
in the previous sections. A more worrisome difference between the different codes is 
the predicted baryon fraction and its profile wi t hin the cluster. Here modern schemes 
still show differences (e.g. see lEttori et alj|2006l : iKravtsov et alj|2005h . which makes it 
difficult to use simulations to calibrate the systematics in the cosmological test based 
on the cluster baryon fraction. 

To date, the comparisons described in the literature show a satisfactory agreement 
between the two approaches, with residual discrepancies originating from the known 
weaknesses which are specific to each scheme. A further limitation of these compar- 
isons is that, in most cases, the simulations are non-radiative. However, at the current 
state of the art, performing comparisons of simulations including radiative losses is not 
expected to provide robust results. As described in the next section, the first relevant 
process that needs to be added is radiative cooling: however, depending on the square 
of the gas density, cooling increases with r esolution without any indication of conver- 
gence, see for example Fig. 13, taken from lBorgani et al.l (|2006T ). At the next level of 
complexity, star formation and supernova feedback occur in regions which have a size 
many orders of magnitude smaller than the spatial resolution of the cosmological sim- 
ulations. Thus, simulations use phenomenological recipes to describe these processes, 
and any comparison would largely test the agreement between these recipes rather 
than identify the inadequacy of the numerical integration schemes. 
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4 Adding complexity 

In this section, we will give a brief overview of how astrophysical processes, that go 
beyond the description of the gravitational instability and of the hydrodynamical flows 
are usually included in simulation codes. 



4.1 Cooling 

We discuss here how the A(u, p) term is usually added in the first law of thermody- 
namics, described by Eq. 1491 and its consequences. 

In cosmological applications, one is usually interested in structures with virial tem- 
peratures larger than 10 4 K. In standard implementations of the cooling function 
A(u,p), one assumes that the gas is optically thin and in ionisation equilibrium. It 
is also usually assumed that three-body cooling processes are unimportant, so as to 
restrict the treatment to two-body processes. For a plasma with primordial compo- 
sition of H and He, these processes are collisional excitation of Hi and Hell, colli- 
sional ionisation of Hi, He I and Hell, standard recombination of Hn, Hell and He ill, 
dielectric recombination of Hell, and free-free emission (Bremsstrahlung). The colli- 
sional ionisation and recombination rates depend only on temperature. Therefore, in 
the absence of ionising background radiation one can solve the resulting rate equation 
analytically. This leads to a cooling function A{u)/p 2 as illustrated in the top panel 
of Fig. 1121 In the presence of ionising background radiation, the rate equations can 
be solved iteratively. Note that for a typical cosmological radiation background (e.g. 
UV background from quasars, see lHaardt fc Madaul fl996l) . the shape of the cooling 
function can be significan tly altered, especi ally at low densities. For a more detailed 



discussion see for example iKatz et al.l 1 19961 ). Additionally, the presence of metals will 



drastically increase the possible processes by which the gas can cool. As it becomes 
computationally very demanding to calculate the cooling function in this case, one usu- 
ally resorts to a pre-computed, tabulated cooling function. As an example, the bottom 
panel of Fig. 1121 at tempera tures above 10 5 K, shows the tabulated cooling function by 
ISutherland fc Dopital 1 19931 ) for different metallicities of the gas, keeping the ratios of 



the different metal species fixed to solar values. Note that almost all implementations 
solve the above rate equations (and therefore the cooling of the gas) as a "sub time 
step" problem, decoupled from the hydrodynamical treatment. In practice this means 
that one assumes that the density is fixed across the time step. Furthermore, the time 
step of the underlying hydrodynamical simulation are in general, for practical reasons, 
not controlled by or related to the cooling time-scale. The resulting uncertainties in- 
troduced by these approximations have not yet been deeply explored and clearly leave 
room for future investigations. 

For the formation of the first objects in haloes with virial temperatures below 10 4 K, 
the assumption of ionisation equilibrium no longer holds. In this case, one has to follow 
the non-equilibrium reactions, solving the balance equations for the individual levels 
of each species during the cosmologica l evolution. In th e absence of metals, the main 
coolants are H2 and hJ molecules ( see I Abel et al.|[l997l ). HD molecules can also play a 
significant role. When metals are present, many more reactions are available and some 
of these can contribute significantly to the cooling function below 10 4 K. This effect 
is clearly visible in the bottom panel of F ig. 1121 for T < 10 4 K. For more details see 



Galli fc PaUal |l99Sl ) or lMaio et all (|2007r i and references therein 
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Fig. 12 The top panel shows the total cooling curve (solid line) and its composition from 
different processes for a primordial mixture of H and He. The bottom panel shows how the 
total cooling curve will change as a function of different metallicity, as indicated in the plot 
(in absolute values). The part below 10 4 K also takes into account cooling by molecules (e.g. 
HD and H2) and metal lines. Taken from Mai o et all (120071 '). 
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Fig. 13 The fraction of cooled baryons / c as a function of the mass of the gas particle, for 
4 different clusters at different resolutions is shown. Filled symbols are for the runs including 
kinetic feedback (e.g. winds), the open circles are re-simulations of one of the clusters with 
wind feedback turned off. The asterisk is for one of the clusters run at very high resolution 
using fewer, but 8 times heavier, gas particles than normal, so that the gas p article mass is 
simila r to that of the DM particles in the high-resolution region. Taken from Borgani ct al. 
J2006f> . 



4.2 Star formation and feedback 

Including radiative losses in simulations causes two numerical problems. Firstly, cooling 
is a runaway process and, at the typical densities reached at the centres of galaxy 
clusters, the cooling time becomes significantly shorter than the Hubble time. As a 
consequence, a large fraction of the baryonic component can cool down and condense 
out of the hot phase. Secondly, since cooling is proportional to the square of the gas 
density, its efficiency is quite sensitive to the presence of the first collapsing small halos, 
where cooling takes place, and therefore on numerical resolution. 

To deal with these issues, one has to include in the code a suitable recipe to convert 
the reservoir of cold and dense gas into collisionless stars. Furthermore, this stellar 
component should represent the energy feedback from supernova explosions, which 
ideally would heat the cold gas, so as to counteract the cooling catastrophe. 

As for s tar fo rmation, a relatively simple recipe is that originally introduced by 
iKatz et al. which is often used in cosmological simulations. According to this 



prescription, for a gas particle to be eligible to form stars, it must have a convergent 
flow, 



Vv; < , 



(98) 
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and have density in excess of some threshold value, e.g. 

Pi > 0.1 atoms cm 3 . (99) 
These criteria are complemented by requiring the gas to be Jeans unstable, that is 

- > 1 , (100) 

where is either the SPH smoothing length or the mesh size for Eulerian codes and 
Cj is the local sound speed. This indicates that the individual resolution element gets 
gravitationally unstable. At high redshift, the physical density can easily exceed the 
threshold given in Eq. 1991 even for particles not belonging to virialised halos. Therefore 
one usually applies a further condition on the gas overdensity, 

> 55.7, (101) 



Pmcan 



which restricts star formation to collapsed, virialised regions. Note that the density 
criterion is the most important one. Particles fulfilling it in almost all cases also fulfill 
the other two criteria. 

Once a gas particle is eligible to form stars, its star formation rate can be written 

as 

dp* _ dpi _ c*Pi nno x 
~dT ~~ ~dt-~' 

where c* is a dimensionless star formation rate parameter and t* the characteristic 
timescale for star formation. The value of this timescale is usually taken to be the 
maximum of the dynamical time fd yn = (47rGpi)~ ' 5 and the cooling time t coo i = 
Wj /(dttj/dt) . In principle, to follow star formation, one would like to produce contin- 
uously collisionless star particles. However, for computational and numerical reasons, 
one approximates this process by waiting for a significant fraction of the gas particle 
mass to have formed stars according to the above rate; when this is accomplished, a 
new, collisionless "star" particle is created from the parent star-forming gas particle, 
whose mass is reduced accordingly. This process takes place until the gas particle is 
entirely transformed into stars. In order to avoid spurious numerical effects, which 
arise from the gravitational interaction of particles with widely differing masses, one 
usually restricts the number of star particles (so called generations) spawned by a gas 
particle to be relatively small, typically 2 — 3. Note that it is also common to restrict 
the described star-formation algorithm to only convert a gas particle into a star par- 
ticle, which correspond to the choice of only one generation. In this case star and gas 
particles have always the same mass. 

To get a more continuous distribution of star particle masses, the probability of 
forming a star can be written as 

p = 1 - exp ^-c*^^ (103) 

and a random number is used to decide when to form a star particle. 

According to this scheme of star formation, each star particle can be identified 
with a Simple Stellar Population (SSP), i.e. a coeval population of stars characterised 
by a given assumed initial mass function (IMF). Further, assuming that all stars with 
masses larger than 8 M© will end as type-II supernovae (SNII), one can calculate 
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the total amount of energy (typically 10 51 erg per supernova) that each star particle 
can release to the surrounding gas. Under the approximation that the typical lifetime 
of massive stars which explode as SN II does not exceed the typical time step of the 
simulation, this is done in the so-called "instantaneous recycling approximation", with 
the feedback energy deposited in the surrounding gas in the same step. 

Improvements with respect to this model include an explicit sub-resolution descrip- 
tion of the multi-phase nature of the interstellar medium, which provides the reservoir 
of star formation. Such a sub grid model tries to model the global dynamical behaviour 
of the interstellar medium in which cold, star-forming clouds are embedded in a hot 
medium. 

O ne example is the multi-phase sub- grid model suggested bv lSpringel fc Hernquistl 
(2003), in which a star-forming gas particle has a multi-phase nature, with a cold 
phase describing cold clouds embedded in pressure equilibrium within a hot medium. 
According to this model, star formation takes place in a self-regulated way. Within 
this model (as within other models), part of the feedback energy is channelled back 
into kinetic energy, effectively leading to a quasi self-consistent modelling of galactic 
outflows, driven by the star-forming regions. Once an efficient form of kinetic feedback 
is included, the amount of stars formed in simulations turns out to converge when 
resolution is varied. An example of the predicted ste llar componen t in dif ferent galaxy 
clusters for varying spatial resolutions, taken from iBorgani et al. I |200d ) is shown in 
Fig. El 

A further direction for improvement is provided by a more accurate description of 
stellar evolution, and of the chemical enrichment associated with star forming particles. 
More accurate models require that energy feedback and metals are released not just 
by SNII, but also by SNIa and low and interme diate mass stars, th ereby avoiding the 
instantaneous recycling approximation (e.g., see lBorgani et aill2008l - Chapter 18, this 
volume, for a more detailed discussion of this point). 



4.3 Additional physics 



A number of other physical processes, besides those related to star formation and SN 
feedback, are in general expected to play a role in the evolution of the cosmic baryons 
and, as such, should be added into the treatment of the hydrodynamical equations. 
For instance, magnetic fields or th e effect of a popu lation of relativistic particles as 
described in detail in the review by iDolag et al. 2008 - Chapter 12, this volume. Other 
processes, whose effect has been studied in cosmological simulation s of galaxy clusters 
include thermal conduction (Jubelgas et aljfeoO^; IDolag et al. 2004), radiative transfer 



to describe the propagation of photons in a medium (e.g. llliev et al.l 20061 and references 
therein), gr owth of black holes and the resulting feedback associated with the extracted 
energy (e.g. ISiiacki et al.ll2007l and references therein. A description of these processes 
is far outside the scope of this review and we point the interested reader to the original 
papers cited above. 



5 Connecting simulations to observations 



Both the recent generation of instruments and an even more sophisticated next gen- 
eration of instruments (at various wavelengths) will allow us to study galaxy clusters 
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Fig. 14 Maps for the SZ decrement for a simulated galaxy cluster. The original map extracted 
from the hydrodynamic simulation, and the same map in the simulated observation (t = 
34 hour) which assumes the AMI interferometric response, are shown in the left and right 
panel, respectively. The side of each map corresponds to 16 ar cmin. The colour scale is shown 
at the bottom of each panel. Taken fromjBonaldi ct al. (20 07T) , 



in rich detail. Therefore, when comparing observations with simulations, instrumental 
effects like resolution and noise, as well as more subtle effects within the observational 
processes, have to be taken into account to separate true features from instrumental 
effects and biases. Therefore, the building up of synthetic instruments to "observe" 
simulations becomes more and more important. 

As an example, Fig. Q3] shows the difference between a Sunyaev-Zel'dovich map 
obtained from a simulated galaxy cluster, and how it would be observed with the AMJj 
instrument, assuming the response according to the array configuration of the radio 
dishes and adding the appropriate noise. Here a 34 hour observat ion has been assumed , 
and the CLEAN deconvolution algorithm applied. For details see iBonaldi et afl(|2007h . 



Another example of a synthetic observation by a virtual telescope is shown in 
Fig. 1151 which shows an opt ical image of a s i mulat ed galaxy cluster with several 
strong lensing features. Here, iMeneghetti et all (|2007t ) investigated the capability of 
the planned DUNE mission to measure the properties of gravitational arcs, includ- 
ing instrumental effects as well as the disturbance by light from the cluster galaxies. 
Shapelet decompositions - based on galaxy images retrieved from the GOODS-ACS 
archive - were used to describe the surface brightness distributions of both the clus- 
ter members and the background galaxies. For the cluster members the morphological 
classifications inferred from the semi-analytic modelling, based on the merger tree of 
the underlying cluster simulations were used to assign a spectral energy distribution 
to each. Several sources of noise like photon noi se, sky backgroun d , seei ng, and instru- 
mental noise are included. For more details see IMeneghetti et al.l (|2007l '). 

Thanks to its high temperature (around 10 8 K) and its high density (up to 10 -3 
particles per cm 3 ) the intra-cluster medium is an ideal target to be studied by X- 



http: / /www.mrao.cam.ac.uk/telescope/ami 
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Fig. 15 Composite ugr + riz + izy image of a simulated galaxy cluster including its simulated 
lensing signal, having its imprint in several strong lensing features. To construct the surface 
brightness distributions of both the cluster members and of the background galaxies, shapelet 
decompositions were used. For each cluster galaxy the morphological classifications and spec- 
tral energy distribution inferred from the semi-analytic modelling based on the merger trees 
from the underlying cluster simulation were used to realistically model the optical appearance 
of the cluster. The fi eld of view is 100" X 100 " and an exposure of 1000 s for each band was 
assumed. Taken from Mcncghct tl et al.l d2007fl . 



ray telescopes. Therefore, so far, most effort towards understanding systematic in the 
observational process of galaxy clusters has been spent interpreting X-ray observations. 

For a direct comparison of the simulations with X-ray observations one has to cal- 
culate from the simulated physical quantities (density, temperature, . . .) the observed 
quantities (e.g. surface brightness). The other way - going from observations to phys- 
ical quantities - is much more difficult and would require a number of assumptions. 
Fortunately, the ICM is usually optically thin, so that absorption of photons within 
the ICM does not have to be taken into account. Hence to obtain an X-ray image of 
the modelled cluster, one has to choose a projection direction and integrate over all 
the emission of the elements along the line of sight for each pixel in the image. The 
X-ray emission at each element is usually approximated as the product of the square of 
the density and the cooling function. As X-ray detectors are only sensitive in a certain 
energy range, one must be careful to take the correct range. For very realistic images 
one needs to taken into account also the effects of the X-ray telescopes and detectors 
(e.g. the limited resolution of the X-ray telescope (point spread function) or the energy 
and angle dependent sensitivy (detector response matrix and vignetting). As observed 
images also depend on the distance of the cluster and the (frequency dependent) ab- 
sorption of the X-ray emission in the Galaxy these should also be accounted for. For an 
exact prediction, the emitted spectrum of each element has to be taken into account, 
as well as many of the effects mentioned above are energy dependent. 
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Fig. 16 Simulated photon images in the 0.7 — 2 keV energy band of a simulated galaxy cluster 
using XMAS-2. The images are binned to 3.2 arcsec. They include background, vignetting 
effects, out-of-time events and the telescope optical paths. From top left to bottom right are 
simulations for the MOS1, PN, and MOS2 instruments on board of the XMM-Newton satellite 
and for th e ACIS-S3 instrum ent onboard of the Cha ndra satellite. Kindly provided by Elena 
Rasia, see iRasia et al.l II2006T ) and iRasia et al.l 1120071 ). 



Even more difficult is the comparison of quantities derived from spectral analy- 
ses such as ICM temperature and ICM metallicity, or their projections in profiles and 
maps. A pixel in an observed temperature map, for example, is derived using all the 
photons within the pi xel area, and the de rived spectrum is fitted with a single model 
of a hot plasma fsee iKaastra et all 120081 - Chapter 9, this volume). These photons, 
however, come from different positions along the line of sight, that have different emis- 
sivities, different temperature and different metallicities. So the spectrum is actually 
a composite of many different spectra. Obviously such a multi-temperature spectrum 
cannot be fitted very well by a single temperature model. Sometimes temperature maps 
are produced in an even cruder way using the ratio of two or more energy bands. These 
are called hardness ratio maps. 

To calculate temperatures, temperature profiles or temperature maps are calculated 
using mostly emission weighted temperatures, i.e. summing the temperature of all 
elements along the line of sight, weighted by their emission 



f WTdV 

T = ^ (104) 

/ WAV K ' 

with T being the gas temperature and W a weighting factor. Usually W is proportional 
to the emissivity of each gas element, 
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W = A(T)n J 



(105) 



where A(T) is the cooling curve and n the gas density. 

This simple procedure does not of course take into account the shape of the spectra 
corresponding to gas of different temperatures. It therefore gives only a rough estimate 
of the actual temperature. With numerical simulations it was investigated how accu- 
rate these emission weighted temperatures are, by comparing them with temperatures 
obtained by actually adding spectra, so-called spectroscopic temperatures. It was found 
that the emission-wei ghted temperatures are systematically higher than the spectro- 
scopi c temperatures ( Mathiesen fc Evrardl 200ll; Gardini et al.l |2004| ; iMazzotta et al.l 
120041 ). To overcome this problem. IMazzotta et alJ (|2004l ) suggested an approximation 
to the spectroscopic temperature, the "spectroscopic-like" temperature T s \ 



Jn 2 T a /T 1/2 dV 
J n 2 T a /T 3 / 2 dV 



(106) 



which yields for a — 0.75 a good estimate of the spectroscopic temperature. In addition 
the sim ulated M-T relati on is strongly affected, if the emission-weighted temperature 
is used (Rasi a et al 

The inhomogenous temperature and metal distri bution was a l so fou nd to be re- 
sponsible also for inaccurate metallicity measurements. iRasia et all ( 20071 ) studied with 
numerical simulations, together with the programme X-MAS2, how well the elements 
Fe, O, Mg and Si can be measured in clusters of different temperature. They find that 
Fe and Si are generally measureable with good accuracy, while O and Mg can be consid- 
erably overestimate d in hot clusters. Using si mulations and the pro gramme SPEX (see 
iKaastra et al.l I2OO8I - Chapter 9 this volume) iKapferer et alJ (j2007T ) found that due to 
the metal inhomogeneities the metal mass in clusters is systematically underestimated 
- in extreme cases by up to a factor of three. 

An example of synthetic X-ray observations is shown in Fig. [16] Here XMAS-2 is 
used to produce photon images from a simulated galaxy cluster. Shown are simulations 
for the MOS1, PN, and MOS2 instruments on board the XMM-Newton satellite, and 
for the ACIS-S3 instrument on board the Chandra satellite. They include instrumental 
effects like back ground, vignetting and response fo r the individual instruments. For 



ettects like background, vignetting and response tor tne lr 
more details see iRasia et all (|2006h and lRasia et~ai1 <|2007h . 



6 Outlook 

In the future, the demand on precision in both simulation techniques and captured 
complexity of the physical processes within the simulations will be quite challenging. 
Recent leading simulations are already extremely difficult to analyse due to their enor- 
mous size and complexity, and they will surely continue to grow. In fact, the next 
generation of supercomputers will grow more in the number of accessable CPUs than 
on the speedup of the individual CPUs and this fact will make the analysis of the 
simulations as challenging as performing the simulations itself. To keep a comparable 
level of accuracy, the interpretation of a simulation of the next generation of high pre- 
cision experiments will need to massively involve virtual telescopes as described in the 
previous section. This will increase the need of involving complex analysis pipelines for 
"observing" simulations, and might lead to a new branch of virtual observers in the 
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astrophysics community, similar to the already, new formed branch of computational 
astrophysicists. 
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